Probabilistic Relative Entropy in Homogenization of Fibrous Metal Matrix Composites (MMCs)

The main aim of this work is to deliver uncertainty propagation analysis for the homogenization process of fibrous metal matrix composites (MMCs). The homogenization method applied here is based on the comparison of the deformation energy of the Representative Volume Element (RVE) for the original and for the homogenized material. This part is completed with the use of the Finite Element Method (FEM) plane strain analysis delivered in the ABAQUS system. The probabilistic goal is achieved by using the response function method, where computer recovery with a few FEM tests enables approximations of polynomial bases for the RVE displacements, and further—algebraic determination of all necessary uncertainty measures. Expected values, standard deviations, and relative entropies are derived in the symbolic algebra system MAPLE; a few different entropy models have been also contrasted including the most popular Kullback–Leibler measure. These characteristics are used to discuss the influence of the uncertainty propagation in the MMCs’ effective material tensor components, but may serve in the reliability assessment by quantification of the distance between extreme responses and the corresponding admissible values.


Introduction
As it is well known, metal matrix composites are composite materials having metal, alloy, or intermetallic matrices [1][2][3], which are reinforced with various types of fibers [4,5], particles [6], or whiskers [7]; layered sheets [8] are alternatively manufactured.Recently, some nanoparticle reinforcements have been invented [9], graphene reinforcements have been implemented [10], and hybrid MMCs have been designed [11].From the material composition point of view, the MMCs are classified into titanium, magnesium, nickel, and copper-based metallic matrices, which offer quite different contrasts of the Young modulus of the matrix with respect to its reinforcement; aluminum matrices are also applicable to minimize the mass of such a composite.Typically, a volume fraction of the fibers reaches 30-40% of the composite, and one of the most popular reinforcing fibers is made of silicon carbide (SiC) [12] or carbon [13]; some alternative reinforcement with carbon nanostructures is also considerable.All these composite materials have been extensively studied in the last 20-30 years due to their wide and important applications in automotive and marine engineering, electronics, and aeronautics.
It is well-known that one of the most efficient mathematical and numerical techniques to analyze composite structures is the homogenization method.Interestingly, this method was and still is applicable to MMCs as the manufacturing method [14].This numerical approach has been applied to investigate thermal stresses [15], MMCs reinforced with short fibers [16], and CNT/MMC under plastic deformation [17].Additionally, the effective thermo-elastic characteristics were calculated [18], and an attempt to analyze MMCs with random fiber distribution has also been described [19].Uncertainty analysis and quantification in the area of MMCs were rather scarce and have been focused on acoustic Materials 2023, 16, 6112 2 of 16 emission [20], damage detection and propagation [21], general reliability analysis [22], as well as strength prediction [23].Analysis of photographic evidence of MMC cross-sections shows that the most important geometrical uncertainties needing attention and numerical simulation are the distribution of the reinforcement into the MMCs domain, interface shape, possible decohesion or lack of contact between the reinforcement and the matrix, as well as reinforcement diameter (and/or length).On the other hand, quite naturally, uncertainties are to be taken into account in mechanical, thermal, and electrical material parameters due to the manufacturing and processing of these composites in higher temperatures and difficulties in obtaining homogeneity within the MMC matrix.
Therefore, the main aim of this paper is to provide probabilistic homogenization [24,25] and uncertainty propagation analysis for some metal matrix composites by taking into account material imperfections within the matrix.Two different numerical techniques are engaged for this purpose, namely the iterative generalized higher-order stochastic perturbation method as well as the semi-analytical technique.They are implemented using numerical recovery of the polynomial bases [26] linking effective elasticity tensor components with the uncertainty source thanks to several FEMhomogenization tests delivered in the ABAQUS system.These techniques enable the determination of the basic probabilistic characteristics of the homogenized tensor such as expectations and standard deviations, while uncertainty propagation is studied thanks to different relative entropy models available from probability theory [27].These include the Kullback-Leibler model [28], Bhattacharyya theory [29], Hellinger idea [30], as well as its Jeffreys symmetrization [31,32].Relative entropy models can be efficiently used to quantify the distance between different pairs of random distributions, which can be represented in structural and computational mechanics by (1) effective and original characteristics of the composites, (2) extreme structural responses and their admissible counterparts (limit functions in the reliability assessment), as well as (3) computed structural stresses and deformations and their counterparts following full-scale experiments.The first case is studied in this paper in the context of all the most popular relative entropy mathematical models recalled above.It needs to be emphasized that all probabilistic mathematical tools together with the Weighted Least Squares Method (WLSM) polynomial bases recovery have been implemented in the computer algebra environment MAPLE, where numerical visualization has been completed.The main motivation of this research was to check how the contrast between the fiber and the matrix elastic moduli affects the probabilistic distance between original and homogenized characteristics of the MMC, which was tested on two different composites, namely MoSiO 2 -SiC and Ti-SiC [24,25].Finally, it needs to be mentioned that this computer strategy may be extended towards thermo-mechanical homogenization of different MMCs, corrosion, and aging stochastic modeling, as well as numerical analysis of the interface defects appearing between various types of reinforcements and metallic matrices.

Homogenization Method
Effective material characteristics of various composites (such as fibrous or particulate, for instance) may be determined using an equation of deformation energies of real and homogenized equivalent structures.The corresponding numerical simulation, most frequently with the use of the Finite Element Method, takes place in the RVE of the given composite domain.Some specific Dirichlet boundary conditions imposed throughout all outer surfaces of this RVE represent uniaxial, biaxial, and transverse deformations, while periodicity conditions are fulfilled on the remaining ones; any von Neumann conditions apply.Using classical notation available in the literature, one writes this energy U (and its equation) as [27] where Ω denotes the RVE (cf. Figure 1), C (e f f ) ijkl and C ijkl are the homogenized and the original composite constitutive tensor, ε ij stands for the strain tensor (adjacent to the adopted geometrical equations), and ε ij represents unitary or zero strains depending on a specific C (e f f ) ijkl component to be determined.where Ω denotes the RVE (cf. Figure 1), ( ) eff ijkl C and ijkl C are the homogenized and the original composite constitutive tensor, ij ε stands for the strain tensor (adjacent to the adopted geometrical equations), and ij ε represents unitary or zero strains depending on a specific  They follow the uniaxial, biaxial, as well as transverse kinematic boundary conditions defined on ∂Ω.The left-hand side integral is determined from the FEM tests in the presence of some uncertainty parameter (b) thanks to the following discretization of the displacements ( ( ) ) and strain tensor ( ( ) where P denotes the order of polynomial basis on Ω, and N stands for the total number of degrees of freedom in the RVE.Let us note that stress tensor discretization is unnecessary to determine homogenized tensor in this approach and is omitted, but possible.Then, the fundamental FEM equations system with applied boundary conditions is solved as where K αβ stands for the global stiffness matrix of the RVE, E is the number of the finite elements in its discretization.A summation of the finite element contributions has no algebraic character except for the one adjacent to FEM global stiffness matrix assemblage; this is performed only once for the entire numerical homogenization process.Therefore, the component They follow the uniaxial, biaxial, as well as transverse kinematic boundary conditions defined on ∂Ω.The left-hand side integral is determined from the FEM tests in the presence of some uncertainty parameter (b) thanks to the following discretization of the displacements (u i = u i (x i )) and strain tensor (ε ij = ε ij (x i )) fields in Ω: where P denotes the order of polynomial basis on Ω, and N stands for the total number of degrees of freedom in the RVE.Let us note that stress tensor discretization is unnecessary to determine homogenized tensor in this approach and is omitted, but possible.Then, the fundamental FEM equations system with applied boundary conditions is solved as where K αβ stands for the global stiffness matrix of the RVE, E is the number of the finite elements in its discretization.A summation of the finite element contributions has no algebraic character except for the one adjacent to FEM global stiffness matrix assemblage; this is performed only once for the entire numerical homogenization process.Therefore, the component C 1111 is determined from Equation (1) while assuming ε 11 = 1 (extension with ∆ 1 and all the remaining components are set to 0), C (e f f ) 1122 -for ε 11 = 1 & ε 22 = 1 (biaxial tension using ∆ 1 and ∆ 2 , zeroes elsewhere), and finally C (e f f ) 1212 -with ε 12 = 1 only (shear of the RVE with ∆ 12 ).Finally: C C These relations include Ω as the total area (or volume in 3D) of the RVE, while Ω f and Ω m denote the domains occupied by the fiber and the matrix, correspondingly.

Probabilistic Governing Equations
Bhattacharyya's relative entropy is a subject of further theoretical and numerical analysis.It quantifies the distance between the probability density functions of the original composite constitutive tensor and its homogenized counterpart.It is possible to calculate this relative entropy with the use of the Bhattacharyya theory as [29] (7) for any subset Ω i of Ω.Some referential models invented by Kullback and Leibler, Jeffreys, and Hellinger have been recalled here: One derives them in the case of two non-truncated Gaussian distributions representing the original fourth-order elasticity tensor (p µ C ijkl , σ C ijkl ) and the effective tensor ) [28][29][30][31][32]: Materials 2023, 16, 6112 5 of 16 A very important case, taking into account the results of further numerical simulation, is when the homogenized tensor is a linear transform of a normally distributed matrix Young modulus, i.e., Then, the first two probabilistic moments of the effective tensor components equal to Further, one derives the following algebraic equations for the aforementioned relative entropies where no summation of i,j,k, or l applies: It needs to be mentioned that all of these formulas include a square of the difference between the expectations of the original elasticity tensor components and their homogenized counterparts.This makes the measurement of the relative entropies of the probabilistic distance more sensitive to this difference than, for instance, the First Order Reliability Method (FORM) distance, where linear interrelation is employed.

Numerical Analysis and Discussion
The Finite Element Method discretization of the RVE consisting of nine parallel and uniformly distributed fibers has been completed in the ABAQUS system; some alternative FEM discretizations and studies can be found in [33,34].It consists of 107.838C3D8 hexagonal finite elements (with a 2 × 2 × 2 integration scheme), where the matrix includes 72.090 bricks, and the fiber-35.748finite elements, cf. Figure 2.This RVE model has been subjected in turn to uniaxial extension in the x direction, to biaxial in the x and y directions (according to notation adopted automatically by the ABAQUS system in Figure 2, which corresponds to x 1 and x 2 from Figure 1), and also to shear deformation in the xy plane to remain in the plane strain state each time.Mean values of material characteristics of the two composites under investigation have been collected in Table 1, where the Young moduli of the matrix have been randomized accordingly.Two sets of initial deterministic FEM analyses have been delivered for the following discrete values of these moduli: 72.090 bricks, and the fiber-35.748finite elements, cf. Figure 2.This RVE model has been subjected in turn to uniaxial extension in the x direction, to biaxial in the x and y directions (according to notation adopted automatically by the ABAQUS system in Figure 2, which corresponds to x1 and x2 from Figure 1), and also to shear deformation in the xy plane to remain in the plane strain state each time.Mean values of material characteristics of the two composites under investigation have been collected in Table 1, where the Young moduli of the matrix have been randomized accordingly.Two sets of initial deterministic FEM analyses have been delivered for the following discrete values of these moduli: (1) Em = {360, 370, 380, 390, 400, 410, 420, 430, 440, 450} GPa, and also ( 2   It has been detected with the Weighted Least Squares Method that polynomial bases linking effective elasticity tensor components with the Young modulus of the metal matrix have a linear form.It has been demonstrated for two MMCs having essentially different contrasts of Ef and Em.This linearity follows a relatively large stiffness of the matrix in MMCs compared to the polymer-based fibrous composites, where polynomial bases have apparently curvilinear character.Therefore, the Second Order Second Moment (SOSM) analysis is sufficient, the Monte-Carlo simulation is unnecessary, and the semi-analytical approach returns simple and exact equations for the moments and the entropies of the homogenized tensor.Therefore, the Gaussian distribution of the Young modulus of the matrix implies Gaussian distribution of the effective tensor and this observation justifies the equations relevant to the relative entropy collected in the previous section.This is a very important result, which remarkably simplifies probabilistic homogenization as no higher moments are necessary and, further, relative entropies have all elegant closed-form algebraic formulas.It is seen that the mechanically driven and polynomial response-based approach is the most efficient in the homogenization of random composites.Probabilistic results of numerical simulations have been collected in Figures 3-20 in two columns, where the left column corresponds to the MoSiO2-SiC composite (with smaller contrast), while the right column corresponds to the Ti-SiC composite.They include expected values (Figures 3-8) and coefficients of variation (CoV, Figures 9-14) of the homogenized tensor as well as relative entropies of this tensor in Figures 15-20.These characteristics have been computed for the components C1111, C1212, and C1122-all as the functions of the CoV of the  It has been detected with the Weighted Least Squares Method that polynomial bases linking effective elasticity tensor components with the Young modulus of the metal matrix have a linear form.It has been demonstrated for two MMCs having essentially different contrasts of E f and E m .This linearity follows a relatively large stiffness of the matrix in MMCs compared to the polymer-based fibrous composites, where polynomial bases have apparently curvilinear character.Therefore, the Second Order Second Moment (SOSM) analysis is sufficient, the Monte-Carlo simulation is unnecessary, and the semi-analytical approach returns simple and exact equations for the moments and the entropies of the homogenized tensor.Therefore, the Gaussian distribution of the Young modulus of the matrix implies Gaussian distribution of the effective tensor and this observation justifies the equations relevant to the relative entropy collected in the previous section.This is a very important result, which remarkably simplifies probabilistic homogenization as no higher moments are necessary and, further, relative entropies have all elegant closed-form algebraic formulas.It is seen that the mechanically driven and polynomial response-based approach is the most efficient in the homogenization of random composites.Probabilistic results of numerical simulations have been collected in Figures 3-20 in two columns, where the left column corresponds to the MoSiO 2 -SiC composite (with smaller contrast), while the right column corresponds to the Ti-SiC composite.They include expected values (Figures 3-8) and coefficients of variation (CoV, Figures 9-14) of the homogenized tensor as well as relative entropies of this tensor in Figures 15-20.These characteristics have been computed for the components C 1111 , C 1212, and C 1122 -all as the functions of the CoV of the Young modulus for the matrix α ∈ [0.00, 0.20].The first two moments have been determined using the 2nd, 4th, 6th, 8th, and 10th order perturbation theories and with the semi-analytical technique.Relative entropies computed here contain the models created by Kullback and Leibler, Bhattacharyya, Hellinger, and Jeffreys.[kPa] [kPa] [kPa]    [kPa] [kPa] [kPa]    [kPa] [kPa] [kPa]    [kPa] [kPa] [kPa]   [kPa] [kPa] [kPa]   [kPa] [kPa] [kPa]                                  The first general observation is that all probabilistic methods return exactly the same first two probabilistic moments for any value of the input α.Additionally, all expectations are constant for any input α, while the resulting CoVs of the effective tensor are linearly dependent on input uncertainty for both composites.Further, as is expected after the data collected in Table 1, expectations of the homogenized tensor for the MoSiO2-SiC composite   The first general observation is that all probabilistic methods return exactly the same first two probabilistic moments for any value of the input α.Additionally, all expectations are constant for any input α, while the resulting CoVs of the effective tensor are linearly dependent on input uncertainty for both composites.Further, as is expected after the data collected in Table 1, expectations of the homogenized tensor for the MoSiO2-SiC composite   The first general observation is that all probabilistic methods return exactly the same first two probabilistic moments for any value of the input α.Additionally, all expectations are constant for any input α, while the resulting CoVs of the effective tensor are linearly dependent on input uncertainty for both composites.Further, as is expected after the data collected in Table 1, expectations of the homogenized tensor for the MoSiO2-SiC composite The first general observation is that all probabilistic methods return exactly the same first two probabilistic moments for any value of the input α.Additionally, all expectations are constant for any input α, while the resulting CoVs of the effective tensor are linearly dependent on input uncertainty for both composites.Further, as is expected after the data collected in Table 1, expectations of the homogenized tensor for the MoSiO 2 -SiC composite are larger than the corresponding characteristics of the Ti-SiC composite; additionally, it has been noted that E C 1122 component uncertainty, where the resulting CoV for the Ti-SiC composite prevails.It is apparent that the resulting uncertainty in the homogenized characteristics is a little bit smaller than the input one.A ratio of this input-to-output uncertainty is not larger than two, with an exception for C (e f f ) 1212 in the case of the Ti-SiC composite, where it is about five.Moreover, while contrasting the CoVs plotted in the left and in the right columns, it is clear that larger uncertainty propagation in homogenization due to Young modulus randomness is noticed for the MMC composite, where The final set of graphs shows a distance of E m probability distribution to C (e f f ) ijkl PDFs for both composite materials.It is evident that all these entropies dramatically decrease while increasing the input CoV.This agrees very well with the fact that increasing these two given PDF ranges makes their distance smaller for a constant difference in their expectations.Moreover, it is documented here that various entropy models return quite different intervals of numerical values.The largest results are obtained while applying the Jeffreys model, smaller numbers are returned with the Kullback-Leibler theory, which is justified by their definitions, followed by Bhattacharyya entropy, and, finally, the Hellinger model; somewhat different relations are exceptionally obtained in the case of the very small distance depicted in Figure 15 1111 , where the first composite's entropy is about four times smaller than the Ti-SiC composite.Summarizing the results concerning relative entropy obtained here, larger uncertainty propagation is expected for the metal matrix composite with smaller contrast between the elastic moduli of this fiber-reinforced composite.The resulting uncertainty in the MMCs' effective characteristics starts to decrease when this contrast increases, and this is demonstrated with an increasing distance of original and homogenized characteristics.Further randomization of Poisson ratios, which is of secondary importance in the MMCs area, may bring larger differences between different probabilistic methods due to its nonlinear impact on the homogenized characteristics.
Let us also note that the relative entropies computed above follow the previous results concerning the first two probabilistic moments so that computational cost has been minimized by an application of the stochastic perturbation technique.This cost is closer to the deterministic FEM solution rather than to the full Monte Carlo simulation results, whereas its accuracy remains the same as for statistical simulations.In summary, it is clearly seen from this analysis that relative probabilistic entropy may be efficient in probabilistic sensitivity analysis, where a smaller distance between two distributions relating some input parameter and the output state function of the given composite is equivalent to a larger impact of this input on the resulting composite's response.

Concluding Remarks
(1) It has been documented in this work that uncertainty analysis in the homogenization of the MMCs brings stable basic probabilistic moments of the effective tensor, which have been computed with minimal time and computer power effort.Uncertainty propagation while randomizing matrix elastic modulus reaches its maximum with E m ≈ E f (for composites with a small contrast between the matrix and the fiber) and decreases while increasing the contrast between these material parameters.Extending this study beyond the limits of the MMCs, rather small uncertainty should accompany homogenized characteristics for the polymer-based composites, where the largest uncertainty is usually observed for the polymers themselves.(2) Relative entropies calculated here due to different mathematical models have enabled the study of the distance between the PDF for E m and the effective tensor components, but the Jeffreys model appeared to be the most distant to the rest of the numerical results.Nevertheless, the determination of these entropies makes it possible to discuss the probabilistic sensitivity of C (e f f ) ijkl components of the given MMC as a function of the aforementioned contrast.This contrast has been discovered as the most influential designing parameter in uncertainty propagation for the given class of composites, so the application of the apparatus presented for the optimization of composite constituents seems to be reasonable and promising.It is documented that future computations of all uncertainty measures for the MMCs may be performed in the future with the lower-order stochastic perturbation method as the fastest approach having the same accuracy as semi-analytical methods, as far as the uncertainty in Young moduli would be considered.This method is relatively easily applicable to any FEM system, contrary to the semi-analytical approach, where symbolic calculus plays a decisive role.The main difficulty while discussing relative entropies for uncertainty of homogenized characteristics is a lack of reference values in the literature and the fact that different theories may return separate numerical value ranges.
(3) Further extensions of this model towards uncertainty quantification in the homogenization process for the nonlinear composites or multi-materials structures subjected to stochastic aging would enable for numerical simulation of some technologically important processes, including some failure and/or corrosion.An application of the stochastic kriging technique [35] or polynomial chaos [36] may be a good alternative to the methods presented above when cross-correlations of multiple uncertainty sources would be expected.An application of the XFEM approach [37] could be beneficial for geometrical randomness homogenization (especially through a few geometrical scales), whereas the probability transformation method (PTM) may enable faster determination of probabilistic (relative) entropy [38].
Funding: This research was funded by research project OPUS no.2021/41/B/ST8/02432 entitled "Probabilistic entropy in engineering computations" sponsored by the National Science Center in Cracow, Poland, 2022-2025.
Informed Consent Statement: Not applicable.

Figure 1 .
Figure 1.Schematic view of the fiber-reinforced composite RVE.

.
The first two moments have been determined using the 2nd, 4th, 6th, 8th, and 10th order perturbation theories and with the semianalytical technique.Relative entropies computed here contain the models created by Kullback and Leibler, Bhattacharyya, Hellinger, and Jeffreys.

Figure 8 .
Figure 8. Expected values of C 1122 for Ti-SiC composite.
both cases.Almost the same conclusion holds true in the case of α C (e f f ) ijkl excluding the C (e f f ) . Comparing the results included in both columns, one notices that the E m distribution is close to C (e f f ) 1111 for the first composite with smaller contrast (E m almost equal E f ).In the case of C (e f f ) 1212 , relative entropies are smaller for the second composite, where E m /E f ≈ 4, but the differences are not so huge as in the first case.The component C (e f f ) 1122 exhibits some similarity to the behavior of C (e f f )

Table 1 .
Material characteristics of the MMCs.

Table 1 .
Material characteristics of the MMCs.